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Summary 


Perturbation  techniques  based  on  Lie  transforms  as  suggested  by 
Deprit  were  used  as  the  theoretical  foundation  for  programming  the 
analytical  solution  of  the  Main  Problem  in  Satellite  Theory  (all 
gravitational  harmonics  being  zero  except  J,) .  The  collection  of 
formulas  necessary  and  sufficient  to  construct  an  ephemeris  is  given 
in  the  exposition.  Short  and  long  period  displacements,  as  well  as  the 
secular  terms,  have  been  obtained  up  to  the  third  order  in  J2  as  power 
series  of  the  eccentricity.  'They  result  from  two  successive  completely 
canonical  transformations  which  it  has  been  found  convenient  not  to 
compose  into  a  unique  transformation.  Division  by  the  eccentricity 
appears  nowhere  in  the  developments — neither  explicitly  nor  implicitly. 
The  determination  of  the  constants  of  motion  from  the  initial  conditions 
has  been  given  an  elementary  solution  that  is  both  complete  and  explicit 
without  being  iterative.  The  program  was  developed  by  Rom  from  MAO's 
package  of  subroutines  for  Mechanized  Algebraic  Operations.  Reliability 
tests  have  been  run  in  two  instances:  in-track  errors  for  ANNA  IB  are 
only  20  cm.  after  210  days  in  orbit,  while  for  RELAY  II,  they  are  2.4  m. 
even  after  350  days  in  orbit. 
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Introduction 


By  the  present  communication,  we  announce  that  the  theory  of  an 
artificial  satellite  can  be  produced  explicitly  in  fully  analytical 
form  by  means  of  programs  which  enable  computers  to  process  literal 
expressions.  We  justify  this  effort  by  displaying  in  two  particular 
instances  the  accuracy  yielded  by  the  series  over  very  long  arcs. 

Theories  of  artificial  satellites  can  be  characterized  by  the  basic 
coordinates  they  use  to  map  the  phase  space.  For  no  other  reason  than 
our  personal  liking  for  Delaunay's  elements,  we  have  chosen  to  develop 
the  Main  Problem  as  set  up  by  Brouwer  (1959).  This  selection  does  not 
imply  a  judgment  on  the  relative  merits  of  Delaunay's  variables  versus 
spheroidal  coordinates  (Vinti,  Kislik,  Aksenov),  elliptic  elements  derived 
from  spheroidal  coordinates  (Iszak)  or  secularly  processing  elliptic  ele¬ 
ments  (Sterne,  Garfinkei,  Aksnes)  .  But  we  submit  that  comparison  of 
analytical  theories  will  not  lead  to  definitive  conclusions  less  we  have 
the  capabilities  of  generating  each  of  them  automatically  by  computer 
so  that  we  can  transfer  analytically  the  constants  of  one  theory  into 
those  of  any  other. 

Our  treatment  of  the  Main  Problem  is  original  on  five  points: 

(i)  We  discarded  the  so-called  Von  Zeipel's  method,  which  is  an 
algorithm  devised  by  Poincare  (1893)  to  generate  a  canonical  transformation 
from  a  determining  function  in  mixed  variables  (old  coordinates  and  new 
momenta).  Instead,  we  use  a  formalism  proposed  elsewhere  (Deprit  1969) 
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under  the  name  of  Lie  transforms.  The  advantage  is  that  we  generate 
explicitly  the  canonical  transformations  and  their  inverses  without 
inversions  or  substitutions,  and  that  we  avail  ourselves  of  a  systematic 
procedure  for  transforming  any  state  function  into  the  new  phase  variables. 

(ii)  We  gave  up  Brouwer's  plan  for  a  theory  in  a  closed  form. 

Indeed,  although  we  could  reproduce  by  Lie  transforms  the  first  order 
terms  computed  by  Brouwer  with  Von  Zeipel's  method,  we  found  that,  at  the 
second  order,  the  quadrature  for  W2  prescribed  by  the  Lie  transforms 
bears,  among  others,  on  terms  of  the  type 


(£-f)r  5sin  2g, 

(£-f)r  4sin(f+2g) 

-  4 

,  -4 

(£-f)r  sin  2g, 

(£-f)r  sir.(f-2g) 

„  -3 

-  7 

(£-f)r  sin  2g, 

(£-f)r  sin(f+2g) 

(£-f)r~3sin(f-2g) 

We  tried  repeatedly  to  express  in  closed  form  the  integrals  of  these 
functions  over  the  mean  anomaly.  Actually,  Moses  (1969)  indicated  that 
such  quadratures  might  not  be  expressible  in  closed  form  by  means  of  the 
usual  elementary  functions.  Similar  terms  have  been  encountered  by  Aksnes 
(1966)  in  his  theory  of  an  artificial  satellite.  In  sum,  the  possibility 
of  implementing  Brouwer's  theory  beyond  the  second  order  in  a  closed  form 
by  means  of  rational,  sine  and  cosine  functions  may  now  be  regarded  as  an 
open  question  that  might  likely  be  answered  in  the  negative.  For  this 
reason,  we  decided  to  expand  the  perturbing  function  in  powers  of  the 


eccentrici tv . 


(iii)  Delaunay's  mapping  from  polar  coordinates  to  elliptic 
elements  has  the  zero  eccentricity  as  a  singularity.  It  causes  troubles  in 
the  corrections  to  be  computed  for  the  mean  anomaly  and  the  argument  of 
perigee.  One  way  of  circumventing  them  would  be  to  coordinatize  the  phase 
space  by  means  of  eccentric  elements  as  defined  by  Poincare  or  Hori.  But 
it  would  generate  considerable,  although  not  insuperable,  complications 

in  the  perturbation  algorithm  (see,  for  instance,  Meffroy  1968).  Instead, 
taking  advantage  of  the  systematic  rule  offered  by  Lie  transforms  for 
transposing  state  functions,  we  decided,  on  one  hand,  to  retain  Delaunay's 
elements  as  the  phase  coordinates,  while,  on  the  other  hand,  we  base  the 
ephemeris  of  the  satellite  on  functions  of  these  phase  coordinates  that 
are  exempt  from  singularities  for  zero  eccentricities.  We  have  selected 
the  mean  distance  F  =  £  +  g  to  the  node,  the  eccentric  functions  C  ■  e  cos  g, 
S  *  e  sin  g,  and  the  usual  Delaunay's  elements,  namely  the  longitude  h 
of  the  ascending  node,  the  polar  component  H  of  the  angular  momentum, 
and  the  action  L  =  via. 

(iv)  We  obtain  the  short  and  long  period  terms  through  the  t hir'd 
order  in  J  ,  and  the  secular  terms  through  the  fourth  order  in  J2.  in 
practice,  so  high  an  order  may  seem  unrealistic.  But,  because  the  series 
are  purely  literal,  they  constitute  a  sort  of  archive  document:  within  the 
accuracy  of  twelve  significant  figures,  their  coefficients  have  been 
determined  once  and  for  all.  In  the  majority  of  cases  where  the  second 
order  is  sufficient,  the  user  can  ignore  in  the  series  the  contributions 

of  the  third  order.  But  for  very  long  arc  predictions,  it  is,  of  course, 
imperative  to  take  into  account  the  secular  terms  of  order  three  and  four. 
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In  sum,  the  higher  the  order  of  a  literal  theory,  the  more  users  it  is  likely 
to  serve. 

Brouwer's  theory  consists  essentially  of  two  successive  canonical 
transformations.  We  tried  to  compose  them  into  a  unique  transformation  as 
Brouwer  proposed  it  in  his  original  paper.  But  the  operations  involved 
so  large  a  number  of  terms  in  the  end  products  as  well  as  in  the  inter¬ 
mediate  results  that  we  concluded  it  would  be  more  economic  to  keep  the 
two  mappings  separate. 

(v)  The  Main  Problem  of  an  artificial  satellite  is  treated  here 
explicitly  as  a  problem  of  initial  conditions.  The  constants  of  the  motion 
are  not  left  to  be  determined  by  successive  iterations  (Cain  1962)  or  by 
least  squares;  the  inverse  canonical  transformations  are  used  to  develop 
explicitly  the  series  that  express  the  average  elements  in  terms  of  the 
osculating  elements. 

The  present  communication  is  well  restricted  in  its  intentions.  We 
meant  more  than  reproducing  by  machine  the  original  paper  of  Brouwer  and 
its  extension  by  Kozai.  In  fact  we  reworked  its  underlying  formalism  and 
saw  to  it  that  it  lended  itself  to  a  smooth  automatic  processing  by 
computers.  At  the  same  time,  we  eradicated  the  main  sources  of  difficulties 
that  so  far  have  adversely  affected  the  use  of  several  theories  of  artifi¬ 
cial  satellites,  namely,  the  premature  truncatures  in  J2 ,  the  determination 
of  the  constants  of  the  motion,  and  the  singularity  of  the  zero  eccentricity. 
Having  established  that  the  Main  Problem  of  Brouwer's  theory  can  be  solved 
to  the  greater  convenience  of  the  users,  we  reckon  that  the  completion  of 


Brouwer's  theory  now  is  more  a  matter  of  developmental  effort  than  a 
research  problem. 

The  basic  steps  of  the  paper  are  the  classical  ones.  We  expounded 
on  the  introduction  of  Delaunay's  variables  with  more  details  than  are 
necessary  for  an  expert  (Section  2);  but,  in  the  classroom  and  in  seminars, 
we  experienced  that  the  standard  references  do  not  adequately  meet  with  the 
demands  of  astrodynamics  engineers  in  this  matter.  Sections  3  -  5  cover  the 
elimination  of  short  period  terms,  i.e.,  the  construction  of  a  canonical 
transformation  to  average  over  the  mean  anomaly,  while  Sections  6  and  7 
outline  the  elimination  of  ;he  long  period  terms,  i.e.,  the  averaging 
over  the  argument  of  perigee  by  constructing  a  second  canonical  transforma¬ 
tion.  We  conclude  the  analytical  study  by  indicating  how  the  secular 
equations  and  the  ephemeris  (in  position  and  velocity)  ought  to  be  computed 
so  that  the  singularity  of  zero  eccentricity  can  be  radically  eradicated 
(Sections  8  and  9  ) , 

Finally  we  present  numerical  evidence  as  to  the  reliability  of  the 
present  solution.  As  illustrations,  we  took  the  satellites  ANNA  IB  and 
RELAY  II  because  their  orbits  have  recently  served  to  compare  various 
theories  (Bonavito  et  at  1968).  Positions  and  velocities  computed  from 
the  series  are  compared  with  the  results  of  an  accurate  stepwise  integra¬ 
tion  of  the  Main  Problem.  The  errors  exhibited  are  unusually  small, 
although  they  are  accumulating  over  time  intervals  comparable  with  the 
period  of  rotation  of  the  perigee.  Our  purpose  here  is  only  to  establish 
that  the  analytical  solution  as  we  have  it  does  not  distort  appreciably, 
even  over  long  arcs,  the  problem  it  purports  to  solve,  which  is  solely 
Brouwer's  Main  Problem  of  Satellite  Theory. 


1. 


The  Main  Problem 


In  the  six-dimensional  phase  space  product  of  the  three-dimensional 
Euclidean  space  of  positions  (x,y,z)  by  the  three-dimensional  Euclidean 
space  of  velocities  (X,Y,Z),  consider  the  Hamiltonian  function 

(1) 

where 

r  *  | x2+y2+z2 | ' 

and  (S  is  the  latitude  with  respect  to  the  coordinate  plane  Oxv, 
thus  unambiguously  defined  by  the  trigonometric  relations 

cos  6  ■  (x2+y2) Vr,  sin  3  •  z/r. 

The  system  described  by  the  Hamiltonian  (1)  constitutes  the  Main  Problem 
(MP)  in  the  theory  of  a  close  satellite  for  an  oblate  planet.  In  that 
context,  m  >  0  is  the  constant  of  gravitation  for  the  planet  (dimension: 
length 3/time2) ,  Rg  >  0  is  the  mean  radius  of  the  planet,  whereas 
Jn  4  0  is  a  (dimensionless)  constant  of  oblateness. 

This  is  a  reversible  system  with  three  degrees  of  freedom.  Its 
Hamiltonian  being  conservative,  it  admits  as  a  first  integral 

m  =  constant.  (2) 

Moreover  the  same  Hamiltonian  being  invariant  with  respect  to  the  commu¬ 
tative  group  of  rotations  around  the  position  axis  Oz,  the  MP  possesses 


y  (X2+Y2+Z2)  -  y 


1 '  * 


(3  sin23-l) 


as  a  second  integral  the  function 


H  =  xY  -  yX,  (3) 

which  is  the  component  along  Oz  of  the  angular  momentum  per  unit  of 
mass  with  respect  to  0  for  the  particle  in  motion  relative  to  the 
coordinate  system  Oxyz. 

The  integral  (3)  is  used  to  render  ignorable  one  coordinate  in  the 
Hamiltonian  (1).  Consider  the  function 

S  S(xNI,yN,h,X,Y,Z) 

-  -x^X  cos  tl+Y  sin  h)  -  y  [  (-X  sin  h+Y  cos  h)2  +  Z2]2. 

For  the  sake  of  brevity  introduce  the  auxiliary  function 

I  I(h,X,Y,Z) 

unambiguously  defined  by  the  consistent  trigonometric  relations 

cos  I  =  (-X  sin  h+Y  cos  h)/[(-X  sin  h+Y  cos  h)2  +  Z2]’4, 
sin  I  =  Z/[(-X  sin  h+Y  cos  h)2  +  Z2]2. 

The  equations 

x  =  -3S/3X  =  x^,  cos  h  -  y„  sin  h  cos  I, 

y  =  -3S/3Y  =  x  sin  h  +  v.,  cos  h  cos  I, 

1\ 

z  =  -3S/3Z  =  y  sin  1, 

i\ 

XN  =  -^S/'-x^,  =  X  cos  h  +  Y  sin  h, 

Y^  =  -3S/r>y^,  =  [(-X  sin  h+Y  cos  h)-  +  Z2  ]  2, 

H  =  -3S/3h  =  xY  -  yX 
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define  implicitly  a  completely  canonical  transformation  from  the  state 
variables  (x,y  ,z  ,X,Y  ,Z)  to  the  state  variables  (x^  ,y^  ,h  ,YN ,H) . 

The  first  three  equations  determine  a  coordinate  transformation  from  the 
position  frame  Oxyz  to  any  frame  Ox^y^N  whose  plane  Ox^y^  contains 
the  position  vector.  This  transformation  composes  a  rotation  about  the 
axis  Oz  with  an  amplitude  h  that  carries  the  axis  Ox  onto  the  axis 
Ox^ ,  and  a  rotation  about  the  axis  Ox^  with  the  amplitude  I  that 
carries  the  plane  Oxy  onto  the  plane  Ox^y^.  Accordingly  the  angle  h 
is  called  the  longitude  of  the  ascending  node  Ox^, ,  and  the  angle  I 
the  inclination  of  the  plane  Ox^.y^  over  the  plane  Oxy. 

The  last  three  equations  of  the  transformation  are  inverted  to  obtain 

that 

X  *  X^  cos  h  -  Y^  sin  h  cos  I, 

Y  =  X^,  sin  h  +  Y^  cos  h  cos  I, 

Z  =  Y  sin  I. 

These  formulas  mean  that  the  plane  Ox^,y^  has  been  selected  to  contain 

permanently  the  velocity  vector,  and  that  (X^.Y^)  are  the  components 

of  the  velocity  vector  in  that  plane.  For  this  reason  the  coordinate 
plane  Ox^y^,  is  cal  led  the  instantaneous  orbital  plane ;  we  shall  denote 
by  H(t). 

The  following  identities  are  trivial  consequences  of  the  transformation 


equations : 
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r  =  x2  +  y2  +  z2  =  x2  +  y2  , 

V2  =  X2  +  Y2  +  Z2  =  +  Y2, 

x  •  X  =  xX  +  yY  +  zZ  =  +  y^. 


The  angular  momemtum  G  =  x  *  X  per  unit  of  mass  is  oriented  along 

‘ 

the  coordinate  axis  ON  normal  to  H(t),  its  component  along  that  axis 
being 


consequent ly 


G  = 


Y.,  - 


Vn 


yZ  -  zY  =  G  sin  I  sin  h, 

zX  -  xZ  =  -G  sin  I  cos  h, 

H  =  xY  -  yX  =  G  cos  I . 

Also  in  the  new  variables  the  formulas  defining  the  latitude  become 

i, 

cos  g  =  (xf.+y2 )  7r,  sin  6  =  (yv/r)sin  I. 

As  a  result  of  the  rotation  from  the  fixed  frame  Oxyz  to  the  moving 
frame  Ox^y^N ,  the  Hamiltonian  of  the  MP  turns  out  to  be 

3t  Ejftx^.-.X^.Y^H) 


In  the  list  of  arguments  for  3i ,  where  the  coordinate  h  was  expected, 
we  place  a  dash  to  emphasize  that  it  has  become  an  ignorable  coordinate. 


For  a  satellite  whose  pi anetocentric  distance  r  is  of  the  order 
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of  R  ,  the  part  of  (4)  which  contains  J,  is  of  the  order  of  J2 . 
Henceforth  assuming  that  is  a  small  number,  we  decompose  (4)  into 

the  sum 

3i  -  Jfn  +  jpfx  (5j) 

with 


In  this  way,  the  MP  is  interpreted  as  basically  a  problem  of  two  bodies 
described  by  (5  ) ,  but  perturbed  by  forces  whose  potential  is  . 

Along  this  line  of  approach  our  interest  lies  in  simplifying  as  much 

as  possible  the  principal  component  J(  .  A  first  step  in  that  direction 

0 

consists  in  introducing  the  polar  components  in  the  instantaneous  orbital 
plane  'I ( t )  .  Thus  consider  the  function 

S  S(r ,”,Xjj,Y^)  «  -rfX^cos  P+YNsin  ) 

and  the  completely  canonical  transformation  from  the  state  variables 
(xN,yN,XN,YN)  to  the  state  variables  (r,0,R/)  implicitly  defined  by 
the  equations 


jy  =  -  tS/  =  r  cos 

yN  =  -•S/,YN  =  r  Sin 
R  =  -  'S /  i r  =  X^cos  0  +  Y  sin  •, 

=  -US/ i"  =  r(-X,,sin  "+Y.,cos  •') . 

M  N 
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We  easily  check  that 

G  =  Vn  '  W 
*»  +  VN  ‘  R2  +  7?  °2. 

so  that  the  Hamiltonian  (5)  transforms  into 

31  :  3f(  r , r  ,- ,R,G,H)  =  JfQ  +  (6^ 

with 

»0  =3fo(r,-.-.R.0,-)  =  7  (R2  +  j?  Q2)  -  7  .  (62) 

X  /i  3  \  3  " 

!  Jf  (r,e,-,-,0,H)  =  -o  “7t|(y  -  2  cos?I)  ■  7  sin21  cos  2eJ*  (63) 

We  aim  at  simplifying  further  the  principal  part  (62)  by  rendering 
ignorable  all  state  variables  but  one  action.  This  ultimate  reduction 
of  the  problem  of  two  bodies  is  accomplished  by  Delaunay's  mapping. 
Consider  the  function 

P  2  P(r ,L,G;u)  =  [-(u2/L2)  +  2(u/r)  -  (G2/r2)f5. 

For  fixed  L  and  G  such  that 

L  >  C-  >  0,  (7) 

P  taken  as  a  function  of  r  has  two  distinct  roots,  namely 
rp  ^  rp ( L , G ; u )  =  [L  -  (L2-G2),2]L/U , 
rA  E  rA(L,G;p)  =  [L  +  (L2-G2)'2]L/;i , 

such  that  r^  >  rp  >  0,  moreover,  P  is  real  under  the  condition  that 
rp  i  r  <  r^,  in  which  case  it  takes  the  form 
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P  =  (u/Lr)[(rA-r)(r-rL)]2. 


Under  the  assumption  (7),  consider  the  function 

S  =  S(r,0,L,G;u)  =  GO  +  J  P(r,L,G;u)dr 

rP 

and  the  completely  canonical  transformation  from  the  state  variables 
(r,0,R,G)  to  the  state  variables  (£,g,L,G)  implicitly  defined  by  the 
equations 


R 

0 

£ 


3S / 9 r  *  P(r ,L,G; u) , 

3S/30  -  G, 

3S/3L  -  (u2/L3)  f  dr/P(r,L,G;p), 
rP 

oS/3C  -  0  -  G  f  dr/r2P(r,L,G;p) . 


It  can  be  expressed  in  a  closed  form.  Indeed  define  the  functions 
a  a(L,u)  and  e  e(L,G)  by  the  conditions 


a  >  0 ,  L‘  =  ua , 

1  >  e  >  0, 

In  terms  of  a  and  e,  the  roots  of  P(r)  =  0  become 


G  =  L(l-e2)'2. 


rP  =  a(l-e). 


=  a(l+e) . 


Now  uniformize  the  quadrature  for  £  by  substituting  for  the  variable 
r  an  angle  u  such  that 


r  =  a( 1-e  cos  u)  ; 
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compute  that 


dr  =  ae  sin  u  du 
-  r  ■  2ae  cos2(u/2), 
r  -  Tp  ■  2ae  sin2(u/2)  , 

and  check  that  the  equation  in  Z  becomes 

i  =  u  -  e  sin  u. 

The  quadrature  defining  g  is  uniformized  by  substituting  for  the 
variable  r  an  angle  f  such  that 

1/r  =  ( 1+e  cos  f ) /p  , 

where  p  p  ( CT ; u )  is  defined  by  the  conditions 
p  •>  0,  G‘  =  y p . 

Calculate  that 


dr  =  (e/p) r ■  sin  f  •  df , 
rA  -  r  =  [ 2e/( 1-e) ] r  cos2(f/2), 
r  -  rp  =  [2e/(l+e)]r  sin2(f/2), 

and  so  check  that 


g  =  •  -  f. 

From  the  resulting  identities 


R  =  (I.e/r)sin  u  =  (ye/f.)sin  f. 
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we  easily  derive  tliat 

R2  +  (C2/r2)  ■=  (2u/r)  -  (u2/L2). 

In  sum,  Delaunay's  mapping  transforms  the  Hamiltonian  (6)  into  the 

sum 


&  -»U,g,-,L,G,H) 


<V 


3f0(-,-,-,L,-,-)  =  -u2/ 2L2, 


(72) 


1  *Re 


=  ^U.g.-.L.G.H)  *  j  r- 


(j  ~  2  cos2lj  -  sin2I  sin(2f+2g)]  (7^ 


2.  The  Case  of  Small  Eccentricities 


We  propose  to  expand  the  perturbation  31  as  given  by  (7  )  in  power 

1 

series  of  the  eccentricity,  thus  limiting  the  application  of  the  theory  to 
close  satellites  with  small  eccentricities. 


The  development  of  1/r3,  cos  2f/r3  and  sin  2f/r3  is  implemented 
automatically  by  computer  (Deprit  and  Rom  1967).  The  following  remarks  are 
in  order.  As  power  series  of  e  with  coefficients  being  trigonometric 
sums  in  the  multiples  of  £,  the  functions  1/r3,  e2cos  2f/r3  and 
e2sin  2 f / r 3  present  the  d'Alembert  characteristic  (Brouwer  and  Clemence 
1961),  Therefore,  if  we  introduce  the  mean  distance 


F  =  l  +  g 


to  the  ascending  node,  we  can  observe  that  the  power  series 


sin(2f+2g) 
_  * 


1  e  j  V 


j-0  k 


A.  ,sin(k?+2F) 
J 
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whose  trigonometric  aiguments  are  not  i  and  g  but  i  and  F  presents 
also  the  d’Alembert  characteris tic ,  i.e.,  the  summation  index  k  satisfies 
both  conditions: 


! k |  <  j  ,  k  i  j  (mod.  2) . 

We  complete  the  expansion  of  by  observing  that 

_i. 

cos  I.  *  H/G  =  (H/l.)(l-e2) 


thus  making  use  of  the  binomial  series  to  develop  cos  I  in  powers  of 
e.  Eventually  trivial  manipulations  of  power  series  in  e  will 
bring  3t j  in  the  form 


u  R*  L 


V  .J 


e 

j“0  L 


1  .j  .0 


+  3f. 


1  »j  »2 


(H/L)- 


where,  for  any  j  2  0,  the  coefficients  Jl  ,  and  are 

1  •  J  » 0  1  i  J  » 2 

finite  sums  of  cosine  functions  with  arguments  of  the  type  pi  +  2kF, 
the  multiples  k  and  p  satisfying  the  conditions 


k  *  0  or  1 ,  [pi  1  j,  p  :  j  (mod.  2). 


To  illustrate  the  development,  we  have  edited  in  Table  I  the  results 
up  to  degree  5  in  e  obtained  by  computer.  We  have  omitted  the  factor 

v. '* R '  L  1  =  n*  a-  (R  /a)"  -  n-  R^. 
e  e 

The  various  trigonometric  arguments  are  entered  in  Column  1;  the  cosines 
are  multiplied  already  by  the  smallest  power  of  e  compatible  with  the 
d'Alembert  characteristic.  Then  the  second  column  lists  the  principal 
parts  of  the  coefficients,  the  next  column  the  parts  of  degree  e’,  and 


so  on . 


r 


Table  1.  Expansion  of  the  Perturbation  Potential 


As  we  shall  see  later,  the  theory  of  satellites  depends  vitally 
on  the  fact  that  the  expansion  of  does  not  contain  the  argument 
2g  *  2F  -  2£.  This  is  proved  to  be  the  case  by  computing  that 

<  *,  >*  ■  ’  1  u'‘^l'V3(1-3  £)• 

The  d'Alembert  character  we  just  emphasized  will  also  prove  quite 
useful.  Indeed  the  perturbation  algorithm  we  are  going  to  use  will  consist 
exclusively  of  Poisson  brackets,  formal  quadratures  and  averagings 
involving  series  that  will  all  derive  from  3f  .  But  the  algebra  of 
series  having  the  d'Alembert  characteristic  is  closed  for  these 
operations.  Thus  by  checking  that  the  initial  input  has  this 

character,  we  make  certain  that  the  for:  ,\lism  will  produce  only  functions 
having  the  same  character.  In  particular,  we  do  not  have  to  make  provision 
for  negative  powers  of  the  eccentricity. 

3.  Elimination  of  the  Short  Period  Terms 

We  propose  to  build  a  completely  canonical  transformation  from  the 
osculating  Delaunay's  state  coordinates  (£ ,g,h,L,G,H)  to  new  coordinates 
(£. '  ,g'  ,h'  ,L'  ,G'  ,H')  such  that,  in  the  transformed  Hamiltonian,  the 
variable  £'  becomes  ignorable.  We  want  the  transformation  in  power 
series  of  J?.  To  this  effect  we  shall  use  the  perturbation  technique 
based  on  Lie  transforms  (Deprit  1969) 

As  we  propose  to  truncate  the  development  after  J^,  we  need  only 


to  fill  in  the  triangle 
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*? 

a"  a!  a§ 

ca/O  <mr  1  m  2  *1/^ 

Jl3  ^t2  1  «-"0 

<uO  nr  1  at  2  —  .  3  fl#  4 

J14  ^«3  Jl2  ™  1 


The  last  line  is  to  be  computed  in  anticipation  of  the  long  period 
elimination  that  we  shall  execute  later  on. 


We  enter  the  triangle  by  putting 


jf°  -  0  for  n  i  2. 

n 


Let  us  indicate  by 


+ 


+ 


the  Poisson  bracket  of  the  functions  <p  and  y  in  the  phase  space 
(£ ' ,g'  ,h ' ,L' ,G' ,M ' ) .  The  calculations  of  the  Poisson  brackets  require 
in  the  present  case  some  attention.  We  have  to  deal  with  functions  whose 
list  of  arguments  are  formally 

L' ,e ' ,H '  '  ,g ' . 

Thus,  on  one  hand,  because  the  variable  h'  is  ignorable  in  both  .;>  and 
1^,  the  calculation  of  (iJv)  will  reduce  to  the  first  two  Jacobians.  On 
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the  other  hand,  L'  appears  In  two  places,  among  the  explicit  variables, 
but  also  implicitly  in  e',  whereas  G'  appears  not  among  the  explicit 
arguments  but  only  implicitly  in  e'.  Therefore  by  differentiation  in 
chain  we  have  that 


36  „  34_  111  .  (l-e'2)^ 

DG 1  “  9e '  3G'  L'e'  ’ 


C8j) 


3L' 


+  Jl_  .  ill 

3e '  SL ' 


(&).■ 


+ 


1  -  e'2 
L'e' 


(82) 


The  square  root  appearing  in  3 / 3G '  will  of  course  be  replaced  by  its 
binomial  expansion  in  powers  of  e*.  Assume  that  the  functions  6 

and  \p  have  the  d'Alembert  characteristic  we  already  mentioned;  then 
the  above  formulas  indicate  that  their  partial  derivatives  with  respect 
to  L'  and  G'  lose  this  characteristic,  and  worse  yet,  contain  terms 
in  e'  *.  Nevertheless,  the  Poisson  bracket  (<{>;ij>)  retains  the  d'Alembert 
characteristic  (Brown  and  Shook  1933).  For  in  the  course  of  calculating 
the  Jacobians 


iL_ilL_il_.i±-  and  il_.lJlL.il_.ilL 

3 £'  DL'  3L '  34'  3g'  3C'  3G'  3g' 

the  terms  in  e'  1  cancel  one  another.  If  the  programming  techniques 
for  manipulating  literal  expressions  by  computers  enable  one  to  process 
polynomial  variables  with  negative  exponents — such  is  the  case  for  MAO 
(Rom  1969) — it  is  relatively  easy  to  write  subroutines  that  will  implement 
the  partial  differentiations  3/9G',  3 / 3L '  as  written  in  ( 8 1 )  and  (8„) . 
Then  the  natural  cancellation  in  the  Poisson  bracket  (fs  ;ij')  of  the  terms 
in  e'  1  may  serve  as  a  check  on  the  validity  of  the  coding  while  in 
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the  initial  phase  of  debugging  the  program,  and  as  a  test  of  the  absolute 
accuracy  on  the  coefficients  obtained  as  soon  as  the  program  becomes 
operational . 

Notice  that  for  any  function  <t>(£ '  ,g'  ,h '  ,L'  ,G'  ,H')  , 

After  these  preliminaries,  we  outline  the  operations  that  accomplished 
the  elimination  of  the  short  period  terms. 

Orde_r  1_.  The  basic  identity  being 

+  o*jj  jWj)  -aj, 

we  selected 

*0  “  i '  * 

Then  putting 

/O  s  <1/0  _  <1/1 
•*1  j  Jtf)  t 

we  obtained  the  generator  Wj  by  the  quadrature 

£  ' 

W,  =  u‘2l  ’3[  3>jdi\ 

Averaging  and  quadrature  are  two  operations  simple  to  program  in  MAO's 
language.  The  average  of  a  function  $  periodic  in  *■ '  is  obtained 
by  transferring  from  the  Poisson  series  s'  into  the  Poisson  series 
<  y  the  terms  whose  trigonometric  argument  does  not  contain  '  ’  .  Then 
the  terms  left  in  constitute  the  properly  periodic  part  of  v;  let 
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us  call  it  J\  As  T  contains  no  constant  term,  the  quadrature 

W  T&l' 

is  a  Poisson  series  having  the  same  type  as  T.  A  term  like  cos(p£'+a) 
of  T  is  transferred  into  W  as  the  term  (l/p)sin(p£'+ct) ;  similarly 
a  term  like  sin(p£'+a)  of  T  is  transferred  into  W  as  the  term 
(-1/p) cos(p£ '+a) .  Table  II  lists  the  development  of  W  up  to  e5. 

The  final  result  is  a  first  order  generator  of  the  form 

W j  -  WjU'.g'.-.L’.C’.H’jR^u) 

=  u2R2L,_3  e,J  [W,  ,  +W  (H'/L')2] 


j^O 


l.j.o  l.j,? 


where,  for  any  j  >  0,  the  coefficients  W  and  W  are  finite 

1 >J  »0  1 »J 

sums  of  sine  functions  in  the  arguments  p£'  +  2kF'  with  k  ■  0  or  1 
and  the  d'Alembert  characteristics  |pj  ^  j,  p  j  (mod.  2).  As  for 
the  first  order  component  in  the  averaged  Hamiltonian,  we  obtain 


where,  for  any  j  -  0,  the  coefficients  J{ 
purely  numerical. 


1 

o  ,j  .0 


and  31 


0 ,  j  ,2 


are 


Incidentally,  until  we  reached  the  end  products  of  the  transformation, 
we  decided  to  keep  explicitly  among  the  polynomial  symbols  of  our  Poisson 
series  the  quantities  p  and  Re>  thus  not  availing  ourselves  of  a 
natural  system  of  units  like  the  Vanguard  units  to  render  dimensionless 
Delaunay's  action  momenta.  Indeed  it  proved  useful  to  constantly  monitor 
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Table  II.  The  first  order  generator  W] 


(The  expansion  must  be  multiplied  by  u2R2l.'  ) 


sin  2F ' 


e'sin  £' 


-  2  +  2  n»2 

8  8  n 


e'sin(i  '+2F' )  |-|  +  |n'2 


e'sin(4'-2F')  _  2  +  2  n'2 


e'-  sin  2t. ' 


8  8  ' 


9  27 

16  ‘  16  n 


e'-sin(2i'+2F')  .11 +  !!„•? 


32  32 


e'2sin(24'-2F') 


e,3sin  3£  * 


MISSING 


53  21  »o 

96  ‘  32  n  “ 


e'3sin(3£'+2F')  ■"4t  +  ^ffri'2 


e’  "sin( 36 '-2F')  ^-^n'2 


e,4sin  kl' 


e'4sin(4£'+2F’)  + 


e*2 

15 

9  ,? 

16 

16  n 

27 

153  l9 

32 

32  n 

123 

-  67  n'2 

64 

64  0 

3  + 

21  n'2 

64 

64  n 

7 

16 

3n'2 

6505  3801  ,, 

1024  ‘  1024  n  ‘ 


.  -2JL  +  _1Z_  ,,.2 

1024  1024 


261  2007  i; 

256  "  256  n 


e'usin(4;.'-2F') 


e 1 1  sin  54' 


e '  4sin(5i ' -2F' ) 


1773 

5319 

2560  ' 

2560  n 

32Gxx  JtU.J.  ,  ; 
5120  5120  n 

81  + 

81 

'  5120 

5120  n 
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the  validity  of  our  codings  as  well  as  the  operational  condition  of  the 
equipment  by  inspecting  the  physical  dimensions  of  the  state  variables 
as  they  were  being  generated.  For  instance,  in  the  present  problem,  all 
elements  of  the  basic  Lie  triangle  have  the  dimensions  of  an  energy 

per  unit  of  mass  (i.e.,  length2 /time2)  whereas  the  generators  have 

those  of  an  action  per  unit  of  mass  (i.e.,  length2 / time) . 

2rieI  !•  Until  the  moment  comes  to  compute  it,  the  second  order 

generator  is  assumed  to  be  the  null  Poisson  series.  Under  this 

assumption  we  calculate  in  the  Lie  triangle  the  provisional  elements 

Si  ♦  «!;“,>  +  «J;w2)  -  «?;»,), 

*0  ‘“I  +  <*»*!>• 

Thus  we  arrive  at  the  second  order  differential  identity 

*1  +  *»o- 

We  select  for  Jf2  the  terms  of  52  independent  from  we  put 

T2  i  ^U'.g'.-.L'.G'.H')  -if2  -3*2 

and  we  obtain  the  second  order  generator  through  the  quadrature 

-?  f1' 

w 2  -  u  L ,37  9  dl\ 

It  turns  out  to  be  a  series  of  the  form 


=  W  (i'.g'.-.L'.G'.H'jR  ,u) 


,W-7Ie'J[W2  +W^  (H’/L’)2+W  u(H'/L')4] 

j>0  2,J,°  2’J’2  2*J»4 


Table  III.  The  second  order  generator  w2 
(The  expansion  must  be  multiplied  by  p 4 R 4 1^ * 


e ' 2 

sin  2F 1 

3  3  ,2  9 

T6  +  8  n  16  n 

459  705  ,2  2457  ,4 

128  32  n  "  128  n 

sin  4F' 

.  JL, 

64  32  1  64 

3  21  ,0  9  ,4 

8  +  32  n  32  n 

e'sin  l* 


e'sin(£'+2F') 


e'sin(H'-2F') 


e'sin(£'+4F') 


e'sin(£ '-4F ' ) 


e'2sin  U' 


e,‘-sin(2£,+2F') 


e'2sin(2£'-2F') 


e'2sin(2£'+4F') 


39  ,  9  ,2  J21  .4 

128  +  64  n  "  128  1 


85  33  ,n  .  179  ,4 

64  "Tn  “  +T4  n 


..135  n, 2  +  225  ,4 

16  n  32  n 


.  .  JLl  n,2  +  15  .4 

132  n  64  n 


9  9  , o  ,  9  .4 

128  "  64  n  +  128  n 


5 _ 63  ,2  _  1227.  ,4 

128  n  256  n 


459  891  ,2  .  1323  ,4 

128  '  ~64  n  +  128  n 


MISSING 


..  IS®  n.2 

128  n 


_1_  +  n'2  . 
128  64  n 


e'2sin(2£.'-4F') 
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where,  for  any  1  £  0,  the  coefficients  W  .  ,  W  and  W  . 

are  finite  sums  of  sine  functions  in  the  arguments  p£'  +  2kF'  with 
k  «  0,  i  or  2  and  the  d'Alembert  characteristic  |p|  -  j*  P  =  j  (mod.  2). 
The  second  order  part  in  the  averaged  Hamiltonian  is  of  the  type 


3f2Q  ■3fJ(-,g,,-,L',G,,H,;Re,p) 


-  li6R  L1"10  Ie'j&f2  .  «£  .  (H'/L')2flg  .  4(H'/L’)4] 

e  j>o  °’j’4 

where,  for  any  j  >  0,  the  coefficients  gfjj  0’&0,j  2  and  ^0  ,j  ,4 

are  sums  of  cosine  functions  in  the  argument  2kg'  with  k  *  0  or  1. 


Order  J3.  The  element  in  the  triangle  can  now  be  completed;  thus 

we  compute 


-#!  -  0>2. 

We  begin  the  calculations  at  order  3  by  putting  W3  ■  0,  and  we  compute 
under  this  assumption  the  elements  in  the  fourth  row  of  the  triangle: 


3t\  “Jf°3  +  +  2C^;WJ  +  ttj; W3)  -  2tH°1;W2), 

Sil  =5f2  +  CtfjjWj)  +  CwJ;W:) 
if03  -aq  +  ewg  jWj) 


Then  we  treat  the  differential  identity 


+  K’V  '*0 

as  we  have  already  done  twice  before.  takes  from  the  terms 

independent  from  2.';  then,  if 

s  _  n/3 
J  ^0  ^0  * 
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Table  IV.  The  third  order  generator  W3 

(The  expansion  has  to  be  multiplied  by  11 


sin  2F ' 

471  ,  8757  JO 741  ,u  .  22455  ,r 

-  1024  +  1024  n  ‘  ‘  1024  n  +  1024  n 

sin  4F' 

sin  6F' 

.  27  +  81  ,2  .  81  „  27  lf) 

1024  1024  1  1024  1  1024  n 

e'sin  £’ 

3  8691  .  12051  ,u  21759  (f, 

52  "  128  n  +  64  n  “  128  1 

e'sin(£'+2F') 

65  ,  10369  „  57677  , .  47243  ,  f. 

512  +  512  n  512  n  +  512  n 

e'sin(£'-2F') 

1089  8973  ,  17163  9279  lfi 

256  -  256  n  +  256  n  *  256  ” 

e'sin(£,+4F') 

1 

■ 

e'sinCi'^F') 

1601  15133  ,o  25463  , 11931  ,r, 

512  512  T|  ,12  '  "  512  ' 

e'sin(£'+6F') 

■  ii  \\m 

e'sin(£’-6F') 

-  27  81  ,  _  81  f4  27  pl6 

512  512  1  512  n  512  1 
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the  third  order  generator  W3  is  derived  from  the  quadrature 

£' 

W3  -  u_2L'3/  ?3di'. 

The  result  is  a  series 


W  =  W  (i'.g'.-.L'.G'.H'ju.R  ) 

>3  j  6 


j6R6L 


.-11 


Ie,J[W  +W  (H'/L')2+W3  .  (H ' /L' ) 1 

j>0  3’j’° 


+W3,j ,6 (»  * /L*  >  6 ] » 


the  coefficients  W,W  ,W  ,W  being  finite  sums  of  sine 

3 » J  *0  3.J.2  3,1.4  3.J.6 

functions  in  the  arguments  p£'  +  2kF'  such  that  0  <  k  <  6,  |p|  <  j 

and  p  =  j  (mod.  2). 


Orde£  As  we  shall  see  later  on,  in  order  to  find  the  third  order 

contributions  to  the  long  period  and  secular  terms,  we  need  to  know  the 
fourth  order  component  in  the  Hamiltonian  averaged  over  the  mean  anomaly 
i.  Consequently  we  have  computed  partially  the  elements  in  the  fifth 
row  of  the  triangle,  namely 

+  3C^;W2)  +30f°:W3)  +  0fO;W4) 

-  30*0  ;W3), 

&22  =3i\  +  CtfiSWj)  +  2(3f};W2)  +  ©fi;W3), 

5f3  * +  Of2;W2); 


then,  arriving  at  the  differential  identity 
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Table  V.  Recapitulation  of  the  Hamiltonian  averaged 
over  the  short  period  angle  and  of  the 
generators  of  the  averaging  transformation. 


W 

1 

*1 

u. 

*5 

W3 

fli  3 

W4 

*0 

e'° 

4 

2 

O 

4. 

6 

3 

12 

4 

20 

5 

e ' 1 

6 

6 

0 

15 

0 

28 

0 

45 

0 

e'2 

8 

6 

2 

18 

6 

36 

8 

10 

e'3 

12 

12 

0 

29 

0 

56 

0 

u 

e'4 

14 

12 

2 

30 

6 

60 

12 

15 

e'5 

18 

18 

0 

45 

0 

84 

0 

0 

e '  6 

20 

18 

2 

45 

b 

84 

12 

20 

e'7 

24 

24 

0 

112 

0 

180 

0 

e'8 

26 

24 

2 

6 

112 

12 

180 

20  | 

e,y 

30 

30 

0 

75 

140 

0 

225 

0  ! 

e»  10 

32 

30 

2 

75 

(> 

140 

12 

225 

20 

e'H 

36 

36 

0 

90 

0 

168 

0 

e.12 

38 

36 

2 

90 

6 

168 

12 

e '  1  3 

42 

42 

0 

0 

e '  1 4 

44 

42 

2 

6 

e.15 

48 

48 

0 

e' 16 

50 

48 

2 

Total 

452 

444 

18 

848 

45 

1200 

72 

1400 

90 
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*S  +  «S:V  -“S- 

.  ++* . 
we  determined  3f4  by  passing  to  it  the  terms  of  3i*  that  do  not 

depend  on  S, 1 . 

The  constructions  we  jast  outlined  have  been  implemented  automatically, 
using  MAO  as  algebraic  processor.  The  development  of  JL  in  power  of  e 
had  been  truncated  after  degree  16,  for  no  other  reason  than  to  set  a 
limit.  At  each  order  of  the  elimination,  two  degrees  in  e  are  lost 
through  the  differentiation  with  respect  to  G’  involved  in  the  Poisson 
brackets.  This  way  we  know  up  to  e16,  but  up  to  e14, 

up  to  e12  and  up  to  e10  only.  By  listing  the  number  of  terms 

in  the  various  components  calculated  so  far.  Table  III  purports  to  suggest 
the  size  of  the  programming  chores  and  how  fast  the  elimination  gains  in 
complexity  as  the  order  Increases. 

4.  The  Elements  of  the  Short  Period  Elimination 

The  generators  Wj,  W2  ,  W3  determine  a  completely  canonical  trans¬ 
formation  from  the  osculating  elements  u  *  (£ ,g ,h ,L,G ,H)  to  a  new  set 
of  variables  u'  ( £ ' ,g' ,h '  ,L'  ,G'  ,H') ,  the  equations  being  of  the  form 

u  -  u'  +  J0u^(u')  +  Tj-  J^u^(u')  +  J^Uq (u ' )  +  ... 

Since  the  coordinate  h'  is  ignorable  in  the  generators, 


H  -  H\ 


On  the  other  hand  the  partial  derivatives  of  the  generators  with  respect 
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to  L'  and  G'  do  not  have  the  d'Alembert  characteristic  in  reference 
to  the  pair  (e * 2. * )  ;  as  a  matter  of  fact,  they  will  contain  negative  powers 
of  e' ,  which  is  a  reflection  of  the  fact  that,  for  circular  orbits, 
the  direction  of  the  perigee  is  undetermined.  But  such  divisors  do  not 
appear  in  the  developments  that  express  the  mean  distance  to  the  ascending 
node 


F  -  £  +  R 


and  the  state  functions 

C  ■  e  cos  g,  S  *  e  sin  g 

in  terms  of  the  new  variables  ( i  '  ,g' ,h ' ,L ' ,C ' ,H ' ) .  Therefore  we  propose 
to  calculate  the  ephemeris  of  the  satellite  by  means  of  the  following 
elements:  F,  h,  C,  S,  L  and  H.  For  all  of  the.;,  except  the  lost  one,  we 
have  obtained  their  expansions  in  powers  of  J  ,  the  coefficients  coming 
in  power  series  of  e'. 

We  review  here  the  formal  characteristics  of  these  series;  the  sta¬ 
tistics  of  Table  VI  make  it  plain  that  we  cannot  think  of  reproducing 
them  in  print.  For  the  sake  of  brevity,  we  have  put 

n'  -  H'/L'  ,  a'  -  a'/R  . 

Also  the  upper  limits  mentioned  in  the  signs  of  summation  are  merely 
anecdotic;  these  are  the  degrees  of  the  eccentricity  e'  after  which  \v 
have  truncated  the  expansions. 

a)  The  mean  distance  F  to  the  ascending  node 
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F  =  F'  +  J  a,_2F'  +  4  J2a'-4F'  +  7  JV_6F1, 

2  1  2  2  2  6  2  3 


F' 

1 

F! 


0<j<14 


0< j <12 


F'  = 

3 


» 

l.J 

+F* 

.0  i,j 

9n'2] , 

1  c 

» 

+F  * 

n’2+Ff 

.j,4 

2.J 

.0  2  ,  j 

3 » j 

+F  *  , 

n '  2+f  ' 

.j.4 

>  0  3  >  j 

,2  3 

,  t  4 


The  coefficients  F'  ,  F'  ,  F*  are  finite  sums  of  sine 

l.j.O’  1 » j  *2  *  3 1  j  » 6 

functions  in  the  arguments  pH *  +  2qF'  such  that  in  F'  we  have 

r  »J 

the  conditions  0  <  q  <  r,  |p|  1  j  and  p  =  j  (mod.  2). 
b)  The  longitude  h  of  the  ascending  node 

h  *=  h'  +  n'U2a'"2h;  +  ~  J2a*~V  +  l  J23a’"6h^), 
h! 


V  p*jh' 

1  4-  ni  1  n» 

0<j<16  1,J»° 


ll2 


I 


0<j<14‘ 


''Xj.o+S.j,/'4' 


h3  *  \  e,j(h^  +h;  .n'-'-Hi;  n'4). 

0<  j<12  •'’J’4 

The  coefficients  h'  ,  h'  .  ,  h'  are  finite  sums  of  sine 

l.j.O  2.J.0  3 ,  j  ,  h 

functions  in  the  arguments  pH*  +  2qF'  such  that  in  h'  ,  we  have 

r  *  j  »k 

0  !  q  i  r,  JpJ  .  j  and  p  j  (mod.  2). 


c)  The  scalar  C  =  e  cos  g 

C  =*  C'  +  J1a'~2C*  +  ^  j£a'”  C'  +  7  J^a,_6C* 

-  1  L  *-  2  O'-  3 


c! 


^  «'j fr' 


0<j<14 


e,J(C'  +c!  .  0n,?) 1 
l.j.O 
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c;  ■  0;|10e'J<c^j.o+c3.J^'’+c’.J^n''+c'-J-6n,r')' 

The  coefficients  Cj^.  C’iJ>2,  ....  C*  J>6  are  finite  sums  of  cosine 
functions  in  the  arguments  pi'  +  (2q+l)F'  such  that,  in  Cr>j>jc  we 
have  Ofqir,  |p|  <  j,  P  =  J  (mod-  2> • 

d)  The  scalar  S  *  e  sin  g 


S  -  S'  +  J0a'”2S'1  +  j  j2at'“4S2  +  ^  S'(  , 

S'  -  y  e' j(S'  +S'  n'2)  , 

1  0<1<14  1*J’°  1,J’2 

s5‘oJu  e'J^.j.o+s;.J,2"'2+s;,-i.."">- 

0<j<10  ,J’  J 

The  coefficients  S ;  _  j  >0  .  SJ  ^  . are  finite  sue*  of  sine 

functions  in  the  arguments  pi'  +  (2q+l)F'  sucii  that ,  in  sr>j^  we 
have  0  i  q  1  r,  |p|  i  j  *  P  j  (rood.  2). 


e)  The  action  L 


L  -  L'(1+J,ci'”JL;  ♦  i  +  j  jy-\), 

v'  '  o;Ii6e,J<l''’J’0+L'’J,2n'2)' 
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The  coefficients  L'  .  L*  L'  .  are  finite  sums  of  cosine 

l.j.O  l.j.2  3 » J  . G 

functions  in  the  arguments  pi'  +  2qF'  such  that,  in  L'  ,  we  have 

r » J 

0  ^  q  1  r,  |p |  f  j  and  p  :  j  (mod.  2). 

Given  the  values  of  F',  h',  C',  S',  L'  and  H'  at  an  instant  t, 
the  preceding  series  enable  us  to  evaluate  the  corresponding  functions 
F,  h,  C,  S,  L  and  H  at  that  instant  of  time.  Indeed  the  decomposition 

pi'  +  2q F'  =  (2q+p) F'  -  pg' 

implies  that  in  all  the  coefficients,  the  elements  e'  and  g'  enter 
only  through  polynomials  in  C'  and  S'.  For  instance,  a  term  like 
e  '  5cos(3i  '+4F' )  can  be  evaluated  as  follows: 

e  '  5cos(  3f.  '+4F')  *  e'-[(e'3cos  3g')cos  7F'  +  (e'3sin  3g')sin  71' 

wherein 

e'-3  -  C'?  +  S'"  , 
e'3cos  3g '  =  4C'3  -  3C' , 
e,3sin  3g'  =  S'(4C':-1). 

Thus  nowhere  is  it  needed  to  evaluate  separately  the  eccentricity  e' 
and  the  argument  of  perigee  g'.  In  sum,  as  far  as  the  short  period 
perturbations  are  concerned,  the  above  series  overcome  entirely  the  indete 
mination  in  the  direction  of  the  perigee  which  is  inherent  to  circular  and 


almost  circular  orbits. 


ri 
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Table  VI.  Recapitulation  of  the  series  necessary  and  sufficient 
to  evaluate  the  short  period  perturbations 
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5.  Initialization  of  the  Elements  F1,  h1,  C*,  S',  L1 

After  the  short  period  terms  have  been  eliminated,  the  Hamiltonian 
becomes  the  function 

M'  -  3V  (~ » g '  »h '  , L '  ,fl '  ,H ' ) 

-  a;  +  -v*;  +  |  Jpe2  *  j  ±  jja-.  (9) 

with 

3i\  3t\  (- »- 1- »L' =  -  7  • 

3i\  =3f;(-,-,-,L,,CMl') 

3^  -»£(-. g'.-.L'.C’.H’)  =  ,  2  <  k  <  4. 

Now  that  is  ignorable  in  the  transformed  variable,  the  action  L' 

turns  out  to  be  an  integral. 

On  one  hand,  through  the  elimination  of  short  period  terms,  the  solu¬ 
tion  of  tii e  Main  Problem  reduces  to  integrating  the  Hamiltonian  equations 
generated  by  (9)  and  composing  their  solutions  with  the  equations  of  the 
canonical  normalization.  On  the  other  hand,  the  initial  conditions  determine 
the  initial  values  of  the  elements  r,  b ,  c,  S,  L  and  H,  but  not  of  the 
transformed  elements.  Therefore,  to  insure  that  the  solution  in 
( i ' , g '  , h ' , L ' , G '  ,  H  ' )  corresponds  through  the  transformation  to  the  initial 
elements  (F^  ,li  ,Cfj  ,S  ,l.,ltH)  ,  we  must  provide  the  way  of  determining,  the 
initial  values  of  the  transformed  elements  (F  '  ,h  '  ,C'  ,S  ' ,  1. '  ,H)  .  In  a 
Calculus  of  Perturbations  based  on  Lie  transforms,  the  initialization  of 
the  transformed  variables  is  accomplished  by  constructing  explicitly  the 
inverse  mapping,  or  rather  the  generators  fV|,V  ,V';)  of  the  inverse 
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mapping.  As  it  was  established  elsewhere  (Deprit  1969),  at  least  up 
to  the  third  order,  they  result  from  these  straightforward  operations: 

Vj  V  (*,*,-tI.,GtH)  =  -W)  (2  ,g,-,L,C,H)  , 

\\  -  V  (5 ,g,-,L,G,H)  =  -W  (£,g,-,L,C,H) , 

‘  t.  J 

V3  e  V3(£,g,-,L,G,H)  =  -W<  -  (Wr  ;Wj)  . 

As  we  did  in  the  case  of  the  direct  transformation,  we  by-passed 
the  explicit  construction  of  the  inverse  transformation  itself,  for  the 
equations  expressing  £'  and  g'  in  terms  of  c,  g,  h,  L,  C,  H  will 
involve  negative  powers  of  e.  Instead  we  addressed  ourselves  directly 
to  the  task  of  expressing  the  elements  F',  h',  C',  S',  L*  in  terms  of 
F,  h,  C,  S,  L.  These  expansions  have  the  same  form  as  the  ones  described 
in  the  previous  section  for  the  direct  transformation.  Inserting  in  these 
series  the  initial  values  FQ ,  hg  ,  Crj ,  SQ,  Lg  and  the  value  for  the 
integral  H  would  then  provide  the  initial  values  Fq  ,  hg  ,  Cq  ,  Sg  and  the 
value  of  the  integral  L'. 


6 .  Elimination  of  the  Long  Period  Terms 


Following  the  main  steps  of  Satellite  Theory  as  outlined  by  Brouwer , 
we  now  consider  constructing  a  completely  canonical  transformation  from 
the  variables  (£  '  ,g' ,h  '  ,L' ,G' ,H  ' )  into  the  variables  (  "  ,g"  ,h" ,1  " ,G"  ,11") 
to  the  effect  of  rendering  the  angle  g"  ignornble  in  the  transformed 
Hamiltonian. 


The  construction  is  again  based  on  determining  the  transformation 
generators  1 , ,  1'  ,  1  . ,  and  so  on.  Thus  entering  a  triangle  similar  to 
the  one  displayed  at  the  beginning  of  Section  3,  we  put 
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3fO(£,,,g",-,L,,,C',,H“)  =^(£”,g’,,-,L,,,G",H"). 
Order  _1.  The  differential  identity 

+  Of0  ;0j)  -aj 

is  satisfied  by  putting 

*2  -a, 

and  by  requesting  that 


3^/3*"  =  0, 

or  that  !  do  not  depend  on  £". 

1 

According  to  what  we  observed  in  Section  3, 


3fj(-,-,-,L",G",H") 


=  u4R‘-L"“6  I  •  (H"/L")r  ] . 

u  »j  »o  o  ,j 


jdO 


Consequent ly 


V/*  _  o  o  • 

2‘"Ren"'-"  1+1  >e"'  W0  ,J+1 ,0^  ,  j+1 /l'")  1 

J-'1 


=  2ii'*RrL""  (l-e":)'*  ^  e"‘ JtW 

P  ' 

j-0 


Mfl 

°tj+l  «0  ()  »J+i 


(H"/L")r ) . 


The  program  has  been  set  up  to  expand  (1-e")2,  to  multiply  this  binomial 
development  with  the  explicitly  written  power  series  in  e"-  ,  and  to 
bring  '&£  /  '>G"  in  a  power  series  of  the  form 
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Da",  C",H") 


3/4u4R2L"‘7A 

e 


1  + 


V  A  e"J 
j>l  J  . 


where 


A 


1 


_  vrl 

5  l..2  » 


and,  for  any  j  -  1,  the  coefficient  is  a  polynomial  in  A  and 

H"/L".  By  introducing  A  as  a  symbol  in  the  list  of  polynomial  variables, 
it  becomes  feasible  to  expand  1/D  in  power  series  of  e"  and,  as  a  matter 
of  fact,  this  artifice  makes  this  inversion  quite  an  elementary  operation 
for  an  automatic  processor  of  Poisson  series  like  MAO.  The  power  series 
1/D  plays  an  essential  rcle  at  the  subsequent  orders  in  the  elimination 
of  the  long  period  term. 


Order  2.  The  element 

3i\  -a®  +  01°;^)  +  C#S;V 

reduces  to 

Jfj 

if  we  assume  that 

3<t2/3£"  -  0, 

i.e.,  if  we  select  for  the  second  order  generator  a  a  function  that 

€. 

does  not  depend  explicitly  on  Consequently 

#2  -3f|  +  ‘3$,  +  2«5;4'1). 

This  relation  constitutes  a  differential  identity  in  the  unknowns 

and  <fj.  We  can  satisfy  by  selecting  for  the  average  of  over 


the  angle  g",  so  that 
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&2Q  e^(-,-,-,L,,,G,,,Hh)  «<#>>„, 


by  assuming  that 


3^/ah"  -  0, 


so  that  the  identity  reduces  to  the  quadrature 


*  2  ’ 


The  results  of  these  operations  are,  on  one  hand  for  the  transformed 
Hamiltonian,  a  second  order  component  of  the  type 

ml  -  mGRV'"10  Ie,,2jpf2  W  .  (H’VL")4] 

u  c  j>q  UfJiO  U  U|J,4 

and,  on  the  other  hand  for  the  canonical  transformation,  a  first  order 
generator  of  the  type 


■I ,  ■  p2R2L"  3sin 


2g"  I  i"je"23  I  «.  .  k 

1>1  Ofkij+1  ’-1, 


(H"/L")2k. 


Ordeir  3.  Entering  the  fourth  row  of  the  transformation  triangle  for 
3C  ,  we  meet  the  element 

m\  +  2 ! £ 2 }  +  Olg S4*3) - 


Assuming  that 


d<t>3/M"  -  0, 


we  reduce  it  to  the  sum 


Jfj  x )  +  2(5^;^). 
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We  isolate  in  it  the  part  yet  unknown  by  putting 

3f\  =J*°3  +  Of2;*x) 

so  that 

3i\  =5^  + 

The  partial  element  is  to  be  computed  whereas  the  Poisson  bracket 

is  set  aside  until  <!>,,  has  been  determined. 

We  treat  in  the  same  way  the  elements 

31]  -#1  +  ; f  1 )  +  CflJ;*.,) . 

»0  +  °?!V- 

By  putting 

-3fJ  +  <*};*,). 

*2  + 

so  that 

3i\  +  3«°:0. 

*0  -50  + 

we  isolate  in  them  the  parts  that  we  are  able  to  compute  and  isolate  the 
Poisson  bracket.  In  fact,  the  relation 


constitutes  a  differential  identity  to  be  satisfied  by  the  unknowns 
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#0?  and 

As  we  have  done  repeatedly,  we  propose  to  take  for  the  average 

of  °ver  g".  and  while  we  assume  that 

3<t„/3h"  -  0, 


we  obtain 


from  the  quadrature 


3$ 


3g"  'HW' 


These  operations  produce,  on  one  hand  for  the  averaged  Hamiltonian,  the 
third  order  component  which  is  a  series  of  the  form 


*o3 


,?R6L 

e 


n-l<.  V 
jiO 


:  a  u 

0sk<1+3  0,:,, 


(H"/L") 


2k 


and,  on  the  other  hand  for  the  canonical  transformation,  the  second  order 
generator 


j‘tR4LM-7V 


jUl 


V  * 

“4+-J  2  » J 


(H,,/L")2k, 


0<k<j+3 


the  coefficients  o  ,  being  finite  sums  of  sine  functions  in  the  argu- 

2  >  J  •  K 

ments  2pg"  with  1  <  p  <  Min(/J,2). 


OnJej;  4^.  Before  computing  the  fourtli  order  row  in  the  triangle,  we 
complete  the  elements  of  the  third  order: 


*i  **?  - 
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Then  in  the  element 

-3*J  +  tt°3 I$x)  +  3CH°;*2)  +  30^;«3)  +  «g;o4). 

we  assume  that 


3$4/3a"  =  o, 


we  isolate  the  terms  that  we  can  already  compute  at  this  stage,  namely 


+  &f3;^  )  +  3 Of2  ;*2) 


so  that 


-3?  +  30^  ;iO 


Similarly  in  the  element 


x\  m3§  +  Cx*;^)  +  2Cxj;$2>  +  C*j;*3> 


we  group  the  terms  that  we  are  already  able  to  evaluate: 


+  CX^j)  +  2«j;*2) 


and  this  way  we  come  to  the  relation 


ari;  -af*  +  40*0^3). 


Proceeding  likewise  for  the  element 


arj  +  «5^2), 


we  put 
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+  #2;^)  +  &f2;4>2) 

and  arrive  at  the  decomposition 

Jfj  =  +  4Gf°;<3). 

Finally,  by  introducing 

5^  =3fJ  + 

we  are  led  to  the  differential  identity 


We  enlist  in  Table  VII  the  number  of  terms  contained  in  the  result 
series  obtained  so  far.  This  summary  indicates  how  the  elimination  of 
short  period  terms  has  simplified  to  a  very  large  extent  the  further 
reduction  of  the  Main  Problem.  But  the  reduction  has  been  bought  at  the 
expense  of  more  sophisticated  manipulations  in  the  construction  of  the 
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r 


now  completely  normalized  Hamiltonian. 

Table  VII.  Recapitulation  of  the  Hamiltonian  averaged 

over  the  long  period  angle  and  of  the  generators 
of  the  averaging  transformation 


*0 

*1 

□ 

■*\ 

»0 

- 

*o4 

e"° 

2 

-  i 

0 

3 

0 

4 

0 

5 

e"2 

2 

3 

3 

5 

5 

7 

7 

e"4 

2 

4 

3 

12 

ft 

16 

8 

e-6 

2 

5 

3 

14 

7 

27 

9 

e"8 

2 

6 

3 

16 

8 

30 

10 

e"10 

2 

7 

3 

18 

9 

33 

11 

e"  12 

*• 

8 

3 

20 

10 

i 

e"14 

2 

9 

3 

e"  16 

2 

Total 

18 

42 

24 

85 

— 

49 

113 

50 

7.  The  Elements  of  the  Long  Period  Elimination  and  Their 
Initialization 

The  generators  of  the  transformation  from  ( i '  ,g' ,h 1 ,1/  ,G * ,H ' ) 
to  ( I" » g" *h",L" ,G"  ,H")  do  not  depend  on  the  angles  and  h"; 

hence  the  transformation  does  not  change  the  actions  I.'  and  H '  : 
thus 


mine 


1/  =  L"  and  H'  =  H”. 

In  order  to  avoid  the  presence  of  divisors  in  e", 
the  transformation  through  its  equations  expressing 


we  do  not  deter- 


h' 


i 
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and  G'  in  terms  of  the  new  variables,  but  rather  by  means  of  the 
elements  FT  *  S, '  +  g'  ,  h',  C'  *  e'cos  g'  and  S’  *  e'sin  g'.  The 
algorithm  by  which  these  state  variables  are  expressed  in  terms  of  the  new 

variables  with  the  help  of  the  generators  0  ,  <f  ,  #  has  been  described 

1  2  3 

elsewhere  (Deprit  1969) .  Here  it  suffices  to  say  to  record  the  morphology 
of  the  resulting  series.  For  the  sake  of  brevity,  we  have  put 


n"  =  h"/l". 


a"  *  a"/R 


a)  The  mean  distance  F'  to  the  ascending  node 


F’  = 


F1 


F"  +  J?a"~2F’1'  +  y  j2a"'4F^  +  j|a"*6F”, 

V  e"2jA_j_1  V  F"  .n"2k, 
l<j<6  0<k<j+2  ,J,k 


V  e"2VJ-2  5 

l-T-5 


f"  n  **  2k 
*2  1  k^  * 

0<k<j+4 


F" 


I  e"2JA-J-3  V  F“  k„ 

1<  j <4  0<k<j+6  ,J,K 


n2k 


The  coefficients  F"  .  ,  are  sums  of  sine  functions  in  the  arguments 

r  •  J  >  k 

2pg"  with  lipi  Min(r, j) . 


b)  The  longitude  h*  of  the  ascending  node 


h'  = 


hl 


h2 


h"  +  n"(J,a"~2hy+  j  J?a"~4h"+  j  jfa,,"6h">  , 

^  i  Z  ^  2  o  ^  3 


V  e,,2jA"j-1  V 


lAj<7 


OikAj+1 

V  e"2jA_J~2  V  h"  .  n"2k, 

l<j<6  0ik<j+3  2,J,k 


U  »» 

Vj.U'1 


2k 
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h"  -  S  e"2j&_j'3  V  h"  .  n,,2k. 

0<T<5  Oikij+5  3,J,k 

The  coefficients  h"‘  .  .  are  sums  of  sine  functions  in  the  arguments 

r,j  ,k 

2pg"  with  1  £  p  <  Min(r,j). 


c)  The  scalar  C'  »  e'cos  g1 


Ci  _  nii  ,  ii  /  ,  ti“2  „  1  t2^ii—  A 

■  C  +  e  (J?a  Cj+  j  JJa  C0+  jr  Jja  , 


2  6  J2l 


CV  -  V  e,,2jA_j“1  I  C".  .  n”2k, 


0<ji6 


Oikij+2 


C"  -  V  e,,2jA"j‘2  y  C'0’  ,  n"Pk, 

0<j<5  0<k<j+4  ’ 

c"  *  V  e"2^A~-^-3  ^  c"  n"2k 

L3  2.  e  a  —  S  4  kn 

0<j<4  0^k<j+5 


The  coefficients  C"  ,  are  sums  of  cosine  functions  in  the  arguments 

r  *j  »k 

(2p+l)g"  with  1  i  p  i  Min(r,j). 


d)  The  scalar  S'  *  e'sin  g* 


S'  -  S"  +  e"(J2a"'2S''  +  \  j2a""4S'2'  +  j  j|ot"_6S”>  , 


S" 


c" 

b3 


y  e.r-)k-i-2  I  S"  , 


S"  -  s  e"23*'3'1  I  *•;  .  kn"2k, 

0ij<6  0ik<j+2 

,..2k 


0£j<5 


0<k<j+4 


V  e"2jA"j_3  V  S"  ,  . n“2k. 


0<j<4 


0<k<j+6 


3 ,  j  »k 


The  coefficients  S"  ,  are  sums  of  sine  functions  in  the  arguments 

f  t  J  »k 

(2p+l)g"  with  1  5.  p  1  Min(r ,  j ) . 


The  inverse  transformation,  i.e.. 


the  mapping  from  (£" ,g" ,h" ,L" ,G" ,H") 


Table  VIII 


Summary  of  the  long  period  corrections 
to  F' ,  h*,  e'cos  g'  and  e'sin  g' 


e” 

AF" 

Ah" 

AC" 

AS" 

0 

0 

0 

3 

3 

2 

4 

3 

8 

8 

4 

5 

4 

10 

10 

6 

6 

5 

12 

12 

8 

7 

6 

14 

14 

10 

8 

7 

16 

16 

12 

9 

8 

18 

18 

14 

_ 

_9 

39 

42 

81 

81 

0 

0 

0 

5 

5 

2 

6 

5 

12 

12 

4 

14 

12 

21 

21 

6 

16 

14 

24 

24  i 

8 

18 

16 

27 

27  j 

10 

20 

18 

30 

30  | 

12 

_ 

20 

74 

85 

119 

119 

0 

0 

0 

7 

7 

2 

8 

7 

16 

16 

4 

18 

16 

27 

27 

6 

30 

27 

40 

40 

8 

33 

30 

44 

44 

10 

33 

89 

113 

134 

134 

Total 

202 

240 

334 

334 
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to  (£' ,g' ,h' ,L' ,C' ,H*)  is  generated  from  functions  1^,  \p2,  i|>  ,  ..., 
which  are  easily  derived  from  the  generators  <?j ,  4>2 »  <t>3 »  •  •  •  of  the 
direct  transformation: 

j  (—  ,g  *  ,L'  ,G'  ,H ' )  , 

^(-.g'.-.L’.C’.H')  *  -^(-.g’.-.L’.G’.H’), 

«3(-.g'. -.L’.G’.H')  -  -*3  -  (<i>2;V. 

These  generators  will  be  used  to  express  the  state  variables  F",  h",  C", 
S"  in  terms  of  F' ,  h ' ,  C' ,  S ' ,  L'  and  H'.  In  particular,  if  we  give 
to  the  elements  the  initial  values  Fg ,  h^  ,  Cg  ,  Sj ,  L'  that  we  computed 

in  Section  5,  we  obtain  from  the  series  the  corresponding  initial  values 
F'g',  hg' ,  Cg  and  S^'  of  the  elements  resulting  from  the  elimination  of 
the  long  period  terms. 

8.  The  Secular  Terms 

After  elimination  of  the  long  period  terms,  the  Hamiltonian  of  the 
Main  Problem  is  the  series 


Jf" 


-3T(-,-,-,L",G",H") 


with 


2  L"2 


ar;  ^;(-,-,-,l",g",h")  -a}. 


a’’  ^:(-,-,-,l",g",h")  -aig. 


The  canonical  equations  it  generates  have  a  trivial  solution.  Since  the 
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angle  coordinates  are  all  ignorable,  the  actions  L" ,  G  and  H  are 
constants  of  the  motion.  Therefore,  the  frequencies 


v 


1 


SL 

3L 


are  also  constants.  Let  us  review  their  morphology: 


a)  Anomalistic  mean  motion  v 


*1  "  n"(1+J2a"  \fl+  \  J2a"  4v!,2+  6  J2“"  6v1,3)» 
Vl.l  *  0<j<7e" 

Le,,2J<Vi.2.j.0+V1.2iJ.2n,,2+V1.2.Ji4f‘",,)' 


1,2 


1.3 


0<j<6 


V  e"2VJ_1  I  v  ,  .  n 

0<k<j+4  1,3,;5,K 


n2k 


0<j<5 


b)  Mean  motion  of  the  perigee  v„ 


c) 


v„  ■  J„n"a"  2(v„  ,+  4-  J_a"  2v.  0+  t-  Jia"  2v„  ,)  , 


-  V  e,,2j(v0  , 

2,1 

0<J<7  2  ’ 1 

-  V  e,,2j(v 

2,2 

o<T<e  2’2 

=  v  e..aJA-J-i 

2  t  3 

0<J<5 

2,1  2  2U  v2 ,2  6  2 

+v>  ,  .  „n”2). 


3 


+v. 


»,  2. 

,n  +v. 


ii4  \ 

n  ) , 


1 1  ■  k 


4  Vo  ,  i  k  1 

0^k< j+4  2,3,J,K 

Mean  motion  of  the  ascending  node  v 


_  .  ,i  m— 2  ti ,  _1  .  m~2 

V  ^  *  J?n  a  h  VV  +  2  v 


3,2  +?  J2“,,‘4v3,3)* 


v3  1  “  1  v3  1  ]e"' ^  • 

'  *  0<3<S  69  1  *J 

—  V  1* 2  j  (  ,  »»2  \ 

V3,2  "  Q<“  ?e  (V3,2,j,0  V3,2,j,2n  * 
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■  1  e  -  -  3  1  kn  * 

Instead  of  integrating  the  differential  equations  in  £."  and  g", 
which  raises  for  the  initial  conditions  i"  and  g'^'  the  problem  of 
divisions  by  the  small  eccentricities  en  and  ej,  we  determine  the 
elements  F",  C" ,  S".  Thus  tlie  integration  of  the  secular  part  of  the 
Main  Problem  results  in  the  following  formulas: 

r  .  +  g-  -  (v1+v:)(t-t0)  +  f”, 

c"  *  e"cos  g"  =  C" cos  v0(t-tQ)  -  S"sin  v2(t-tQ), 

S"  *  e"sin  g"  *  S"cos  v,(t-t0)  +  C"sin  v2(t"t0)» 
h11  =  v,(t-tc)  +  h|\ 

9 .  Computation  of  the  Coordinates  and  Velocities 

After  the  series  for  the  direct  and  inverse  transformations  eliminating 
all  periodic  terms  have  been  obtained,  we  are  in  a  position  to  compute  at 
any  instant — within  the  time  interval  of  validity  of  the  theory — the  position 
and  the  velocity  of  the  satellite  from  coordinates  and  velocity  components 
at  an  epoch  t  . 

Given  (x  ,y,,z  ,  x  ,v  ,  z(1)  »  the  initial  values  of  the  elements  F() , 
hp ,  C>,  S  ,  L., ,  H  are  determined  from  the  following  system  of  formulas: 

The  components  of  the  angular  momentum  being 


A'  *  yoz-i  -  zny  , 


B  =  z-Xo  -  x-jZr,, 

H  ■  Vo  -  Vo- 


we  compute  its  norm  at  t  , 
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An  +  Bn  +  H 


so  that  the  trigonometric  relations 


Ao  *  Go  sin  lo  sin  h0* 


Bq  ■  -G0  sin  I0  cos  h0. 


H  *  G0  cos  lo 


determine  unambiguously  the  initial  longitude  hc  of  the  ascending  node 
and  the  initial  inclination.  Next  we  rotate  the  coordinate  axes  so  that 
O'  and  On  lie  in  the  orbital  plane  n(tQ)  with  0£  coinciding  with  the 
ascending  node  at  the  initial  instant  t  .  In  this  frame  of  reference, 
the  initial  coordinates  and  velocities  of  the  satellite  are 


«o  ■  xo  cos  ho  +  sin  V 

n0  -  (-xQsin  h0+yQcos  h0)cos  IQ  +  zQ  sin  IQ , 

*  *0  cos  ho  +  sin  ho  * 
nn  *  (-x0sin  h0+y0cos  h0)cos  I0  +  z0  sin  I0. 

Evidently  the  initial  planetocentric  distance  is 


r0  "  ^0  +  V 


Moreover,  from  I.aplace-Hamilton's  vector  oriented  along  the  apsidal 
line,  we  obtain  that 


co  "  co^o  “  forn 


i  =  “co'-o  “  r,cro  * 


which  relations  determine  the  square  of  the  eccentricity 


eo  *  cb  +  s‘o’ 


Having  evaluated  the  quantity 
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2E  =  (SJ+ng)  -  2r^, 


we  are  now  in  a  position  to  calculate  the  initial  semi-major  axis 


a  -  -1/2E  , 

0  0  * 


hence  Delaunay's  action  at  time  t  ,  namely 


L0  =  '/V 


There  remains  to  compute  F  .  Recalling  that 


—  =  COS  -  C„  + 

a 


^0  "  co  +  - - 7==  (c0sin  *cTV°s  ♦o>* 

1  +  Vl-et 


—  =  sin  6.  -  S„  - 

a  y0  o 


1  + 


—==  (CQsin  *0-S0cos  *„>, 


where  =*  EQ  +  gQ ,  (En  being  the  eccentric  anomaly  at  time  tQ)  ,  we 
deduce  that 


cosin  *0  •  Socos  *0  cono  -  Vo 


1  +  /l-e2 

so  that  the  preceding  formulas  yield 


a0/l-e2 


cos  *0  =  C0e2  +  —  + 


sin  ^  =  Sce2  +  —  - 


^0  .  Wo“soV  1  “  *  1_eo 


Wo-V0>  *  “  J  1_eo 


y/1^2 


They  determine  unambiguously  ij>0 .  Finally  an  obvious  modification  of 
Kepler's  equation  produces 


Fo  =  ^0  -  c0  sin  *0  +  SQ  cos  V 

The  initial  values  Fr  ,  h(1 ,  CQ ,  SQ,  LQ ,  H  are  entered  into  the 
series  expressing  the  averaged  elements  F' ,  h' ,  C' ,  S ' ,  L'  in  terms  of 
the  osculating  elements  F,  h,  C,  S,  L  and  H;  the  substitution  produces 
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the  initial  values  F’,  hg  ,  Cg  ,  Sg  and  the  value  of  the  integral  L*. 

They  are  then  introduced  in  the  series  expressing  F",  h" ,  C"  and  S" 
in  terms  of  F' ,  h ' ,  C' ,  S 1 ,  L'  and  H,  which  results  in  the  determination 
of  the  initial  values  Fg,  hg,  dg,  and  Sg . 

At  this  stage  we  are  in  a  position  to  evaluate  from  the  series  the 
numerical  values  of  the  three  basic  mean  motions  v  ,  v2  and  v3.  This 
in  the  last  step  in  the  initialization  phase  of  the  ephemeris. 

The  calculation  of  the  position  and  velocity  at  any  instant  t  follows 
the  same  line  as  the  initialization,  but  in  the  reverse  order. 

First  we  evaluate  the  state  variables  F",  h",  C"  and  S"  at  time 
t  from  the  simple  formulas  given  at  the  end  of  Section  8.  Then  the  series 
described  in  Section  7  provide  the  average  state  variables  F’,  h',  C*  , 

S',  L'  at  time  t,  whereas  the  series  of  Section  4  furnish  the  osculating 
elements  F,  h,  C,  S,  L  at  that  instant. 

After  the  process  of  evaluating  the  series  numerically  is  completed, 
we  determine  the  coordinates  and  components  of  the  velocity  from  the 
following  system  of  formulas: 

Kepler’s  eauation  in  the  form 

,  *  F  +  C  sin  .  -  S  cos  , 


is  solved  by  iteration  to  obtain  the 
and  the  semi-major  axis  are  computed 

e'  = 


anomaly  c  *  E  +  g.  The  eccentricity 
from  the  relations 


Cr  +  S: 


a  *  L-  . 
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The  position  in  the  nodal  frame  of  reference  is  given  by 


E,  »  a^cos  ip  -  C  + 
n  ■  a|sin  ^  -  S  - 
Knowing 


(C  sin  \p-S  cos 


*>] 


1  +  V  1-e2 

C  (C  sin  ip— S  cos  t^)j  . 


1  +  y/  1-e2 


G  -  L  /l  -  e2  , 


r  -  /c2  +  n2  » 


we  can  use  Hamilton's  vector  to  get  the  components  of  the  velocity  in  the 
nodal  frame, 


?  -  -  +  r)' 

*  *Hc  +  f)- 


Thereafter,  having  produced  the  inclination  from  the  relations 


cos  I  ■  H/G, 


0  <  I  <  ir, 


we  perform  the  usual  rotations  to  pass  from  the  nodal  frame  to  the  original 
inertial  frame,  and  we  come  finally  to  the  Cartesian  coordinates 


x  *  £  cos  h  -  n  cos  I  sin  h, 
y*£sinh  +  n  cos  I  cos  h, 
z  ■  n  sin  I, 


and  the  components  of  the  velocity 

x  *  5  cos  h  -  n  cos  I  sin  h, 
y  *  £  sin  h  +  n  cos  I  cos  h, 

7  n  sin  T 


10.  Reliability  Tests 

While  the  coding  was  in  development,  repeated  checks  were  applied  to 
test  its  correctness.  The  series  generated  by  the  program  were  constantly 
checked  for  their  physical  dimensions  and  their  d'Alembert  characteristic. 
Subroutines  that,  in  the  parlance  of  the  trade,  are  called  expert ,  like 
taking  partial  derivatives  with  respect  to  Delaunay's  actions  G  and  L 
or  computing  a  Poisson  bracket,  have  been  thoroughly  tested.  For  Instance, 
in  regard  to  the  canonical  transformation  (£,  ,g,h,L,G,H)  -*•  (£'  ,g'  ,h'  ,L'  ,G'  ,H')  , 

we  satisfied  ourselves  that  the  Poisson  bracket  (£:g)  is  indeed  equal 

-12 

to  zero,  i.e.,  that  the  coefficients  in  the  result  are  smaller  than  10 
in  relative  accuracy.  (N.B.  They  should  have  come  exactly  equal  to  zero 
had  we  operated  in  Integer  arithmetic.)  The  general  course  of  a  reduction 
by  Lie  transforms  has  been  implemented  on  a  number  of  simple  examples  like 
the  simple  pendulum,  Duffing's  equation  and  the  relativistic  harmonic 
oscillator. 


The  ultimate  test  of  reliability  is  a  comparison  of  the  positions 
and  velocities  predicted  by  the  series  with  those  of  a  highly  accurate 
numerical  Integration  scheme. 

We  chose  to  integrate  the  equations  of  the  Main  Problem  by  recurrent 
power  series.  All  Taylor  expansions  involved  are  computed  at  each  step 
through  degree  16;  the  time  step  is  selected  so  as  to  maintain  12  significant 
figures  in  both  the  integral  of  energy  and  the  integral  of  polar  angular 
momentum.  The  procedure  is  rapid,  highly  accurate,  and  stable;  it  enables 
one  to  follow  for  very  long  arcs  the  orbit  corresponding  to  a  given  set  of 
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initial  conditions  (Deprit  and  Zahar  1966). 

In  order  to  allow  comparisons  with  other  prediction  methods,  we 
borrowed  from  Bonavito  et  al  (1968)  the  initial  Cartesian  coordinates 
and  velocities  they  assign  for  the  artificial  satellites  RELAY  II  and 
ANNA  IB.  Table  IX  lists  the  initial  conditions,  the  osculating  elements 
at  epoch,  the  corrections  to  be  applied  in  order  to  obtain  from  them  the 
constants  of  the  motion  as  per  our  theory,  and  the  basic  periods. 

We  examined  the  residuals  (integration-series)  for  the  elliptic  ele¬ 
ments  selected  in  our  theory.  Disagreements  on  the  inclination  I  and 
in  the  longitude  of  the  node  h  are  simply  insignificant.  As  expected, 
the  mean  distance  F  to  the  node  shows  a  secular  deviation  (see  Fig.  1) . 
For  ANNA  IB,  the  residuals  in  the  semi-major  axis  suggest  a  long  period 
error  coupled  with  a  secular  trend;  for  RELAY  n ,  the  presence  of  a  long 
period  error  is  well  marked  in  the  element  C.  For  all  other  variables 
presented  in  Figure  1,  the  effects  of  short  period  errors  caused  by 
the  truncatures  in  the  eccentricity  and  seem  to  mask  long  period  and 

secular  tendencies. 


More  revealing  quantities  to  evaluate  long  range  reliability  are 
the  intrinsic  deviations.  If  x^,  y^,  z ^  (resp,  yg,  zg)  are  the 
coordinates  at  time  t  furnished  by  the  numerical  integration  (resp. 
the  literal  series),  and  if  X^,  Yj,  Z^  (resp.  Xg,  Yg,  Zg)  are  the  compo¬ 


nents  of  the  velocities,  the  directions  of  the  tangent,  binormal  and 
normal  to  the  orbit  at  time  t  are  given  by  the  triplets 
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Table  IX.  Constants  of  the  motion  for  the 
test  orbits 


Initial  osculating  elements 
FQ(in  radians) 
hQ(in  radians) 

so 

Co 

LQ(in  Vanguard  units) 

H  (in  Vanguard  units) 

Short  period  corrections 


Long  period  corrections 


ANNA  IB 


RELAY  II 


2.538 

875 

214 

278 

3.273 

083 

992 

516 

0.949 

636 

751 

294 

-2.384 

959 

105 

384 

-0.002 

107 

639 

831 

-0.025 

229 

668 

345 

-0.006 

371 

881 

838 

-0.234 

623 

580 

641 

1.085 

131 

662 

111 

1.322 

050 

356 

567 

0.695 

348 

576 

283 

0.884 

318 

864 

870 

0.273 

044 

549 

x  10“ 3 

-0.052 

347 

711 

x  10 

0.342 

375 

395 

x  10-3 

-0.006 

726 

021 

x  10 

-0.369 

163 

708 

x  10“  3 

0.123 

600 

234 

x  10 

0.015 

809 

392 

x  10“3 

0.563 

272 

260 

x  10 

-0.128 

216 

782 

x  10“  3 

-0.452 

874 

015 

x  10 

0.005 

752 

x  10 

-6 

0.846 

017 

x  10 

-6 

0.012 

755 

x  10 

-6 

2.070 

715 

x  10 

-6 

-0.561 

488 

x  10 

-6 

-2.404 

673 

x  10 

-6 

1.410 

317 

x  10 

-6 

21.619 

075 

x  10 

-6 

Periods 
2tt/v  i 
2tt  h2 
2tt/|v3| 


ln47m10?2139 

121d05h38m15?35 

99d17h34m55!60 


3h15"\)0!4528 
331d03h50m54?43 
332d22h12m5 7*11 


Millimeters 


M*on  Dittanc*  to  th*  Mod*  F 


30  «  10-" 

10  >  lO"5' 
0 

-1.0  K  10""  | 

-3.0  «  10-” 


100 

Days 


200 


Eccentric  El*m*nt  C  =  *cotg 


ANNA  1-B 


RELAY  II 


Fig.  1.  Errors  on  osculating  elliptic  elements 
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II  "  ‘WWW 

where 

vi  “  <xi+yi+zi>,>. 
bj.  -  (Aj/GIiBI/GI,CI/GI) 

where 

^  “  yizi  -  ziV 
Bi  "  zixi  ■  xizi* 
ci  *  xiyi  ■  yixi* 

Gj  -  (Aj+Bj+Cj)^, 

and 

nx  -  ((BIZI-CIYiyGIVI,(CIXI-AIZI)/GIVI,(AIYI-BIXI)/GIVIj. 

Then  the  dot  products 

QsrV  •  h>  (ir^s)  '  Zi'  (ir^s}  *  £i 

constitute  the  projections  of  the  error  in  position  respectively  on  the 
tangent  ("in-track  error”)!  the  normal  ("along-track  error")  and  the 
binomial  ("across-track  error")  of  the  orbit  at  time  t.  It  is  a 
characteristic  feature  of  a  perturbation  theory  that,  while  it  yields  a 
very  close  approximation  of  the  orbit  even  over  very  long  range,  it 
leaves  in  error  the  basic  clocks  of  the  motion,  namely  the  frequencies 
v i ,  V  ,  v 3 .  However  small  the  level  of  errors  on  these  clocks,  they 
cumulate  linearly  with  time.  The  predicted  orbit  coincides  nicely  with 
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the  real  trajectory  of  its  initial  conditions,  but  the  satellite  takes  a 
more  and  more  pronounced  habit  of  arriving  too  early  or  too  late  at  the 
expected  meeting  places.  This  story  can  be  read  in  the  diagrams  of  position 
errors  (Figure  2).  The  deviations  along  and  across  the  track  are  totally 
insignificant,  but  the  errors  in  track  show  a  secular  trend.  (Note  that 
for  both  satellites,  the  tests  cover  more  than  one  revolution  of  the  peri¬ 
gee.)  Comparison  with  similar  tests  by  Arsenault  et  al  (1963)  ,  Lubowe 
(1966)  and  Bonavito  et  al  (1969)  should  restore  confidence  in  the  capa¬ 
bilities  of  a  satellite  theory  based  on  Delaunay's  elements.  For  ANNA  IB, 
after  more  than  3000  revolutions,  the  in-track  error  reaches  only  0.2  meter; 
as  for  RELAY  II,  after  2800  revolutions,  it  is  still  as  small  as  2  meters. 

Conclusions 

In  Perturbation  Theory,  Lie  transforms  will  likely  supersede  Von 
Zeipel's  method;  they  provide  easy  routine  schemes  for  inverting  canonical 
transformations,  for  determining  the  constants  of  the  motion,  and  for 
transposing  state  variables. 

Analytical  expansions,  however  large  the  number  of  terms  in  the  series, 
are  capable  of  very  high  accuracy  over  long  intervals  of  time.  As  a  matter 
of  fact,  a  Satellite  Theory  carried  through  the  fourth  order  in  its  secular 
terms  can  deliver  the  accuracy  presently  contemplated  by  radar  and  laser 
experiments. 

Computers  can  be  programmed  to  generate  developments  of  average 


complexity. 


Anna  1-B 


Across -Track 


Relay  II 


Intrinsic  errors  in  position  (integration-series) 
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The  present  theory  eradicates  the  difficulties  caused  by  small 
eccentricities.  It  is  evidently  incomplete.  But  it  has  solved  in  principle 
the  problems  to  be  encountered  from  including  more  terms  of  the  gravi¬ 
tational  field,  the  luni-solar  perturbations  and  other  perturbations  derived 
from  force  functions. 

Readers  interested  in  having  the  series  produced  by  the  present 
algorithm  should  contact  the  first  author  of  the  paper. 
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